Demo 3 - Remote parallel inversion [distributed]

In [1]:

    

# profile = 'phobos'   # remote workstation
# profile = 'pantheon' # remote cluster
# profile = 'zeus' # remote workstation
profile = 'mpi' # local machine


    

In [2]:

    

import numpy as np
from zephyr.Dispatcher import SeisFDFDDispatcher
from zephyr.Parallel import CommonReducer
from IPython.parallel import Reference


    

In [3]:

    

import matplotlib.pyplot as plt
import matplotlib.cm as cm
import matplotlib
%matplotlib inline

import mpld3
mpld3.enable_notebook()


    

In [4]:

    

lclip = 2000
hclip = 3000
clipscale = 0.1
sms = 0.5
rms = 0.5

def plotField(u):
    clip = clipscale*abs(u).max()
    plt.imshow(u.real, cmap=cm.bwr, vmin=-clip, vmax=clip)

def plotModel(v):
    plt.imshow(v.real, cmap=cm.jet, vmin=lclip, vmax=hclip)

def plotGeometry(geom):
    
    srcpos = geom['src'][:,::2]
    recpos = geom['rec'][:,::2]
    
    axistemp = plt.axis()
    plt.plot(srcpos[:,0], srcpos[:,1], 'kx', markersize=sms)
    plt.plot(recpos[:,0], recpos[:,1], 'kv', markersize=rms)
    plt.axis(axistemp)


    

In [5]:

    

cellSize    = 1             # m
nx          = 164           # count
nz          = 264           # count
freqs       = [1e2, 2e2, 3e2, 4e2] # Hz
freeSurf    = [False, False, False, False] # t r b l
nPML        = 32            # number of PML points
nky         = 1            # number of y-directional plane-wave components


    

In [6]:

    

velocity    = 2500          # m/s
vanom       = 100           # m/s
density     = 2700          # units of density
Q           = 500           # can be inf


    

In [7]:

    

srcs        = np.array([np.ones(101)*32, np.zeros(101), np.linspace(32, 232, 101)]).T
recs        = np.array([np.ones(101)*132, np.zeros(101), np.linspace(32, 232, 101)]).T
nsrc        = len(srcs)
nrec        = len(recs)
recmode     = 'fixed'

geom        = {
    'src':  srcs,
    'rec':  recs,
    'mode': 'fixed',
}


    

In [8]:

    

cache       = False         # whether to cache computed wavefields for a given source
cacheDir    = '.'

parFac = 1
chunksPerWorker = 0.5       # NB: parFac * chunksPerWorker = number of source array subsets
ensembleClear = False


    

In [9]:

    

dims        = (nx,nz)       # tuple
rho         = np.fliplr(np.ones(dims) * density)
nfreq       = len(freqs)    # number of frequencies
nsp         = nfreq * nky   # total number of 2D subproblems
cPert       = np.zeros(dims)
cPert[(nx/2)-20:(nx/2)+20,(nz/2)-20:(nz/2)+20] = vanom
c           = np.fliplr(np.ones(dims) * velocity)
cFlat       = c.copy()
c          += np.fliplr(cPert)
cTrue       = c


    

In [10]:

    

# Base configuration for all subproblems
systemConfig = {
    'dx':   cellSize,       # m
    'dz':   cellSize,       # m
    'c':        c.T,        # m/s
    'rho':      rho.T,      # density
    'Q':        Q,          # can be inf
    'nx':       nx,         # count
    'nz':       nz,         # count
    'freeSurf': freeSurf,   # t r b l
    'nPML':     nPML,
    'geom':     geom,
    'cache':    cache,
    'cacheDir': cacheDir,
    'freqs':    freqs,
    'nky':      nky,
    'parFac':   parFac,
    'chunksPerWorker':  chunksPerWorker,
    'profile':  profile,
    'ensembleClear':    ensembleClear,
#     'MPI': False,
#    'Solver':   Reference('SimPEG.SolverWrapD(scipy.sparse.linalg.splu)'),#Solver,
}


    

In [11]:

    

%%time
sp = SeisFDFDDispatcher(systemConfig)
sp.remote.dview.activate()
sp.remote.e0.activate('e0')
survey, problem = sp.spawnInterfaces()
srcs = survey.genSrc()
sp.srcs = srcs


    

    




CPU times: user 244 ms, sys: 33.1 ms, total: 278 ms
Wall time: 855 ms

In [12]:

    

%%time
sp.forward()


    

    




CPU times: user 26.7 ms, sys: 3.16 ms, total: 29.8 ms
Wall time: 30.9 ms

In [13]:

    

sp.forwardGraph


    

    
Out[13]:

In [14]:

    

%%time
dObs = sp.dPred
sp.dObs = dObs


    

    




CPU times: user 3.71 s, sys: 763 ms, total: 4.47 s
Wall time: 30.6 s

In [15]:

    

%%time
sp.rebuildSystem(cFlat.T)
sp.forward()
dPred = sp.dPred


    

    




CPU times: user 3.92 s, sys: 812 ms, total: 4.73 s
Wall time: 32.7 s

In [16]:

    

dResid = sp.residual
misfit = sp.misfit


    

In [17]:

    

%%time
sp.backprop()
g = sp.g


    

    




CPU times: user 3.38 s, sys: 738 ms, total: 4.12 s
Wall time: 29.6 s

In [18]:

    

%pylab inline


    

    




Populating the interactive namespace from numpy and matplotlib

In [19]:

    

ind = dResid[0].real
clip = abs(ind).max()
imshow(ind, vmin=-clip, vmax=clip, cmap=cm.bwr)
colorbar()


    

    
Out[19]:





<matplotlib.colorbar.Colorbar instance at 0x1132f4b90>






    

In [21]:

    

ind = (g[0]+g[1]+g[2]+g[3]).real
clip = abs(ind).max()
imshow(ind, vmin=-clip, vmax=clip, cmap=cm.bwr)
colorbar()


    

    
Out[21]:





<matplotlib.colorbar.Colorbar instance at 0x1138a0ef0>






    

In [ ]: